In [1]:
import numpy as np
import matplotlib.pyplot as plt

In [2]:
plt.rcdefaults()
# Typeface sizes
from matplotlib import rcParams
rcParams['axes.labelsize'] = 12
rcParams['xtick.labelsize'] = 12
rcParams['ytick.labelsize'] = 12
rcParams['legend.fontsize'] = 12
#rcParams['font.family'] = 'serif'
#rcParams['font.serif'] = ['Computer Modern Roman']
#rcParams['text.usetex'] = True

# Optimal figure size
WIDTH = 350.0  # the number latex spits out
FACTOR = 0.90  # the fraction of the width you'd like the figure to occupy
fig_width_pt  = WIDTH * FACTOR

inches_per_pt = 1.0 / 72.27
golden_ratio  = (np.sqrt(5) - 1.0) / 2.0  # because it looks good

fig_width_in  = fig_width_pt * inches_per_pt  # figure width in inches
fig_height_in = fig_width_in * golden_ratio   # figure height in inches
fig_dims      = [fig_width_in, fig_height_in] # fig dims as a list

rcParams['figure.figsize'] = fig_dims

In [3]:
%matplotlib inline

In [ ]:
import pandas as pd
pd.set_option('display.mpl_style', 'default')
import itertools


/Users/jcole/anaconda/lib/python2.7/site-packages/pytz/__init__.py:29: UserWarning: Module argparse was already imported from /Users/jcole/anaconda/lib/python2.7/argparse.pyc, but /Users/jcole/anaconda/lib/python2.7/site-packages is being added to sys.path
  from pkg_resources import resource_stream

First, let's load the results from the small model of polled included in the default settings. This involves loading four animal files (live cows, dead cows, live bulls, and dead bulls). We will load them and merge them into a single data frame.


In [ ]:
for sim in xrange(1,11):
    # Load the individual history files
    fij = pd.read_csv('holstein/1/fij_paa_pryce_%s.txt'%sim, sep=' ', header=None, names=['calf','sire','dam','gen','fij','paa','mating'])
    if sim == 1:
        all_animals = fij
    else:
        all_animals = pd.concat([all_animals, fij])

In [ ]:
all_animals.head()

In [ ]:
fij_on_paa = pd.ols(y=all_animals['fij'], x=all_animals['paa'], intercept=True)
print(all_animals['fij'].corr(all_animals['paa']), fij_on_paa.f_stat['p-value'])